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Abstract. Polaron and bipolaron formation in the Holstein-Hubbard model 
with harmonic confinement potential, relevant to quantum dot structures, is 
investigated in one to three dimensions by means of unbiased quantum Monte 
Carlo simulations. The discrete nature of the lattice and quantum phonon effects 
are fully taken into account. The dependence on phonon frequency, Coulomb 
repulsion, confinement strength (dot size) and electron-phonon interaction 
strength is studied over a wide range of parameter values. Confinement is found 
to reduce the size of (bi-)polarons at a given coupling strength, to reduce the 
critical coupling for small-(bi-)polaron formation, to increase the polaron binding 
energy, and to be more important in lower dimensions. The present method also 
permits to consider models with dispersive phonons, anharmonic confinement, or 
long-range interactions. 
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1. Introduction 

Continuous advances in fabrication methods in semi-conductor physics over recent 
decades have enabled researchers to create nanoscale systems in which carrier motion 
is confined along one or more spatial dimensions. Among these structures are quantum 
dots [1,2], corresponding to quasi-zero-dimensional systems with some atom-like 
properties [3], which are of substantial technical interest due to their potential use 
as lasers [4], in quantum computing [5], as storage devices [6], or as single-photon 
sources [7, 8] . Improvements in experimental techniques now also permit studies of 
the properties of individual dots rather than ensembles [9], as well as of molecular 
quantum dots [10]. 

Apart from strong correlations between electrons, the interaction of charge 
carriers with the lattice degrees of freedom is of great relevance, especially concerning 
its role in dephasing and relaxation processes [11-13] as well as transport [14]. 
Evidence for the existence of polarons (carriers bound to a self-induced lattice 
distortion) has been given [15], and non-trivial multi-phonon effects and polaron 
formation have been proposed to explain the lack of experimental evidence for the 
phonon bottleneck in quantum dots [11, 16]. Formation of bipolarons (bound electron 
pairs sharing a phonon cloud) has been proposed as an explanation for pair tunneling 
from GaAs quantum dots [17]. Besides, (small) bipolaron states significantly influence 
the transport through single molecules [18, 19]. 
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The problem of polaron formation in a quantum dot has received a lot of attention 

in the past (sec [20 -24] and references therein). From this work, the general conclusion 
is that confinement enhances the tendency of an electron to undergo a crossover to a 
(small) polaron so that heavy polaronic quasiparticles may be realised experimentally 
even in the weak or intermediate elcctron-phonon coupling regime. Less work has 
been devoted to understand bipolaron formation in quantum dots [17,23-26], with 
contradictory results on the effect of confinement [23, 26]. 

All existing work is based on the continuum Frohlich model [27] , employing mainly 
variational approaches of different reliability — the above-mentioned conflicting results 
seem to originate from the inadequacy of some of the approaches used. However, from a 
theoretical point of view, the most interesting case is that in which the extension of the 
lattice displacement attached to the charge carriers is comparable to the confinement 
length. The Holstein molecular-crystal model has been originally developed [28] to 
answer the question as to whether a local lattice instability occurs upon increasing 
the electron-phonon coupling. This problem cannot be addressed in the framework of 
continuum models, since the local lattice dynamics on the scale of the unit cell has to 
be taken into account [29]. 

The need for lattice models also stems from the fact that it is nowadays possible 
to fabricate self- assembled quantum dots with lateral dimensions of less than 4 nm [9] , 
which contain only a small number of unit cells in each direction. Obviously, in such 
systems, the discrete nature of the lattice will play an important role, motivating us 
to revisit the problem of (bi-)polaron formation. 

The Holstein-Hubbard model considered here is not completely realistic, as it 
only includes local electron-phonon and electron-electron interactions as well as a 
coupling to dispersionless optical phonons. However, we are interested in fundamental 
effects arising from the combination of electron-phonon interaction and confinement 
in discrete systems. Besides, many aspects of quantum dots may be understood using 
effective models [1]. 

From numerous studies of Holstein models with one (two) electrons on a discrete 
lattice [30,31], it is well known that the non-linear process of (bi-)polaron formation 
represents a complex many-body problem that may not be satisfactorily described by 
means of variational or perturbativc approaches. The necessity for a universal all- 
coupling approach to describe intermediate-coupling materials such as CdTe has been 
pointed out long ago [21], motivating the application of unbiased numerical methods. 

In this work, we extend the world-line quantum Monte Carlo (QMC) method of 
[32-35] to include a harmonic confining potential. The method yields practically exact 
and unbiased results in any dimension, in principle with no restrictions of system size 
or parameters. As a consequence, we will be able to study (bi-)polaron formation over 
the whole range from weak to strong confinement. The method is capable of treating 
systems with long-range interactions or dispersive phonons, which opens up interesting 
perspectives for future work. An important finding as compared to previous work 
(including [36]) concerns the problem of very long autocorrelations times in certain 
parameter regimes. Finally, in contrast to earlier Green-function QMC calculations 
for a continuum model [17], the present approach does not rely on a trial wavefunction. 

The paper is organised as follows. In section 2, we present the model and briefly 
discuss the underlying approximations. The QMC method is described in section 3, 
and section 4 contains a discussion of our results. We end with a summary and an 
outlook to future work in section 5. 
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2. Model 

We consider a Holstein-Hubbard model, supplemented by a harmonic confining 
potential. Using dimensionless phonon variables, the Hamiltonian reads 

^ = E 4s. + c/ E + ^ E + Y E + E • (1) 

(i^j)a i i i i 
^— — V ' ^ ^ ' ^ ^ ' ^ ^ ^ ' 

Hkin ffee ifoon -f^ph=-f^ph + ^ph -^'P 

Here cj^ creates an electron of spin cr on lattice site i, and -ffi^-j^ describes the hopping 
of electrons between nearest-neighbouring lattice sites {i,j) with hopping integral t. 
The second term, H^^, accounts for on-site Coulomb repulsion between electrons, 
where n^^ = c\^c^^. The harmonic confinement potential, its strength being measured 
by is given by H^^^, r- denotes the position vector of an electron at site i, and 
"•i = Scr^icr- ^ph constitutes the kinetic (-ffp^) and potential energy (H^yi) 
lattice degrees of freedom, corresponding to independent harmonic oscillators with 
energy uIq. Finally, H^^ mediates the coupling between the local electron occupation 
n- and the lattice distortion x^ with coupling strength a. We use units such that 
h = = e = 1, consider D-dimensional hypcrcubic lattices of linear size A'' and 
volume N^, and assume periodic boundary conditions in real space. 

In terms of the coupling parameter a, the atomic-limit {t = 0) polaron binding 
energy Ep is given by Ep = 0^/(210^). We define the usual dimensionless coupling 
constant A = 2Ep/W, with the free bandwidth W = 4tD. Furthermore, we introduce 
the adiabaticity ratio 7 = tOf^/t, which permits to distinguish between the adiabatic 
(7 < 1) and the non-adiabatic (7 > 1) regimes. 

This paper is exclusively concerned with the cases of either one electron, or two 
electrons of opposite spin forming a spin singlet. Obviously, for N^ = 1, H^^ = and 
spin indices may be dropped. In experiments, the maximum number of electrons in a 
quantum dot is strongly influenced by the dot size (set here by K), so that the low- 
density regime considered here may be relevant in particular for small dots, which are 
of interest due to their potential use as quantum bits [6] . One and two-electron states 
also play an important role in molecular transistors, corresponding (in the simplest 
case) to a quantum dot with a single electronic level coupled to a vibrational mode 
[18,19]. 

As we are interested in fundamental effects due to confinement, the Holstein- 
Hubbard-type model seems a good starting point, especially as a lot of knowledge is 
available for K = [35,37 40]. Although quantitative results for realistic quantum 
dot systems are beyond the scope of this work, let us briefly discuss the simpliflcations 
inherent to the Hamiltonian (1). 

Using a one-band model we ignore the existence of a band gap in semi-conducting 
host materials. Hence, our considerations are restricted to one or two carriers in 
the same (either conduction or valence) band. A symmetric, parabolic confining 
potential — which can be realised experimentally [1] — is assumed for simplicity, 
although more general situations may be studied. 

We restrict ourselves to a simple effective electron-phonon coupling to a 
dispersionless (optical) phonon branch. For such local phonon modes, and in the 
absence of sharp interfaces, the bulk-phonon approximation is expected to be a good 
starting point. The coupling to optical phonons has been shown to dominate bipolaron 
formation [17]. 
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The Hamiltonian (1) only contains on-site interactions, which are not completely 
justified in the case of semi-conducting materials. Long-range interactions will be the 
subject of future investigations. 



3. Method 



The QMC method used here is a straight-forward extension of [32,33,35]. The 
appealing features of this approach are the analytical integration over the phonon 
degrees of freedom — enabling us to study the adiabatic regime — and the fact that 
the numerical effort is essentially independent of system size. Besides, the formalism 
treats the cases of one and two electrons on the same footing, and is general enough 
so as to permit future studies of models with, e.g., dispersive phonons or long-range 
interactions. 



3.1. Partition function 

The derivation of the relevant fcrmionic contribution (the bosonic part can be 
calculated exactly) to the partition function Z = Tr e~^^ is almost identical to [33, 35]. 
Dividing the imaginary-time axis [0,(3] (/? = (fegT)"-^ is the inverse temperature) into 
intervals of length At = (}/L we may write 

L 



Tr 



g-Ar(/f,<i„-FHco„)e-Ar(H-h-FHep)e-ArHJ^ 



The result for the partition function reads [33] 

with the fermionic weight 
^f({r«)})=exp| Pir-r') £ 



(2) 
(3) 

(4) 



r,r'=l 



X exp • 



r=l 



U „(2) 



L We 
T=1J=1 



£=1 

Here N^ = 1 (2) for the polaron (bipolaron) problem, and the components of the 
position vector of electron ^ on time slice r, r'^^K arc denoted as r^^^j^ (fi = 1, . . . , D). 
The fermion world-lines are subject to periodic boundary conditions both in real space 
and imaginary time, and the sum in equation (3) is over all allowed configurations. 

The retarded electron (self-) interaction due to electron-phonon coupling is 
described by the memory function 



4L 1 - cos(27ri//L) + {u^^Tf 12 

and electron hopping contributes 

J(r) = ^ E '^^^(^ • ^) • 



(6) 
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We define the expectation value of a static observable O as 

(O) = Z-'TrOe-''" = Z^' ^ 0({r(«)})«;f ({r(«)}) . (7) 

The fermion contribution to the total energy is obtained from 

= — ^ In Zi = El +E^ +E^ + E^ (8) 

L kin ' ep ' ee ' con V"/ 

with the kinetic energy {S runs over all nearest-neighbour sites) 
and the interaction energies 



t,t'=i 5,5'=i ^ ' 

^ee = X^E(ia,,.<^)), (11) 
r=l 

and 

EL = jKj:f:{\r^Of). (12) 

Also of interest is the electron-lattice correlation function measured in the 
direction /z = 1, 

N 

Cep(r) = {EpPNj-' ' r = 0,l,...,N-l, (13) 

which fulfills the sum rule Y.r <^ep(^) = 1 (^^^ Ep>0). Within QMC, we have [34] 

L We 

Cep(r) = {E^PK)-' Yl - E('5r-'«>,r'1) ) ' (1^) 

r,T'=l £=1 " 

Similarly, for N^ = 2, we calculate the electron-electron correlation function 

N 

CeeW = E("iTW' r = 0,l,...,N-l, (15) 

i=l 

fulfilling C'ee('') = 1) for which the QMC estimator reads 
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3.3. Simulation details 

The system described by the partition function (3) is characterised by an additional 
dimension (imaginary time), as well as by a comphcatcd retarded interaction. As 
first shown in [32], it may be simulated by means of Markov Chain Monte Carlo in 
combination with the Metropolis-Hastings algorithm [41]. To improve convergence for 
critical parameters, wc use both local updates (change of a random component of the 
position vector on a random time slice of one particle by one lattice site) and global 
updates (translation of an entire world line by one lattice site) [35]. 

An important point overlooked in previous work [32 34,36] is the; problem of 
significant statistical correlations between successive configurations due to the local 
updating. Whereas previous authors performed measurements every L steps, we find 
that this is by far not sufficient to ensure statistically independent measurements, 
especially at intermediate electron-phonon coupling. The integrated autocorrelation 
time in such cases may well exceed tens of thousand of MC steps, so that a careful 
binning analysis [42] is required for every run to ensure correct results and error bars. 
Details on this problem will be reported elsewhere. Although autocorrelation times 
are expected to be generally shorter in continuous-time simulations (see, e.g., [43]), 
the problem must not be neglected. A QMC algorithm entirely free of autocorrelations 
has been presented in [40,44], but this approach suffers from restrictions in system 
size for D > 1 [45]. Another problem of world-line algorithms is that the acceptance 
rate for local updates approaches zero in the strong-coupling regime. For K/t > 1 
and A > 1, this becomes also true for global updates. Finally, we also detected 
convergence problems in the vicinity of the small-(bi-)polaron crossover due to critical 
slowing down. 

The only systematic error in the present calculations comes from the Suzuki- 
Trotter approximation (equation (2)). It may be eliminated by working in continuous 
imaginary time [43], or by performing simulations at different At and exploiting the 
Ar^-dependence of the results to scale to At = 0. To save some computer time, here 
we have simply chosen a single (small) value of At = 0.05, which ensures that the 
Suzuki- Trotter error is satisfactorily small. 

4. Results 

We discuss our numerical results in relation to the available knowledge about the 
Holstein-Hubbard model (equation (1) with K = 0). Similar to the non-confined case, 
the adiabaticity ratio 7 turns out to have an important effect on the physics, and 
a comparison of the cases 7 <C 1 and 7^1 will be given for the ID polaron. In 
higher dimensions, we shall focus on the adiabatic regime 7 ^ 1, as the low- lying 
optical phonon energies in quantum dot materials are in the range of 10-100 meV, i.e., 
substantially smaller than the electron bandwidth. 

The 2D (3D) geometry considered here — corresponding to a circular (spherical) 
quantum dot -may actually be realised experimentally, whereas the ID case is mainly 
of theoretical interest. However, we shall see that the physics is qualitatively similar 
in D = 1 and D > 1. 

All simulations have been carried out using a linear lattice size of = 31, enough 
to obtain very well converged results even for if = [33,36]. Since for if > 
local physics becomes more important, finite-size effects are negligible. The center 
of the parabolic potential is located at site i = of each dimension, and the site 




Figure 1. (a) Kinetic energy Ekin of one electron in a ID system in the adiabatic 
regime (7 = 0.1) as a function of electron-phonon coupling strength A for different 
values of the confinement strength K. Here and in subsequent figures N = 31, 
fSt = 15 and Ar = 0.05. Lines are guides to the eye only, (b) Electron-lattice 
correlation function Cep (r ) as a function of distance r for the same parameters as 
in (a) and A = 1. The inset shows Cep(O) as a function of A. 



indices range over [—15,15]. To study ground-state properties, we set the inverse 
temperature = 15. The number of Trotter time sUces has been fixed to L = 300 
so that At = 0.05. Errorbars arc typically smaller than the linewidth, and are shown 
only if larger than the symbols used. 

A drawback of the present approach (cf [43,46]) is that we can not calculate the 
quasiparticlc effective mass. Instead, wc shall consider the electronic kinetic energy, 
which includes contributions from incoherent processes. Additionally, the correlation 
functions defined in section 3 permit us to monitor the size of the (bi-)polaron, and 
to determine the critical coupling for the crossover to a small-(bi-)poIaron state. 

The confinement length — often a directly tunable parameter in previous work — is 
set by the oscillator strength (~ K) of the harmonic potential. For two electrons, it 
may be estimated from the electron-electron correlation function C^^{r) at A = 0. 

4-1- Polaron 

Existing work suggests that the ground state of the Holstein model with one electron 
{K = 0, U irrelevant) at weak electron-phonon coupling is a large polaron (extending 
over more than one lattice site) in ID, and a quasi- free electron for D > 1 [47]. 

With increasing coupling A, the potential energy due to lattice deformation 
increases relative to the kinetic energy of the electron. In the adiabatic regime " < 1, 
a small-polaron state (with the electron and the lattice distortion being localised 
essentially to the same site) is formed if Ep > W/2 (A > 1). In contrast, for 7 > 1, 
small-polaron formation requires a sizeable lattice distortion, leading to the condition 
g"^ = Ep/ujQ > 1. Note that these "critical couplings" are conceptually different from 
those for the occurrence of the Peierls quantum phase transition at half filling [48] 
since there is no phase transition in the (bi-)polaron problem for 7 > [49]. 

The small-polaron state is characterised by a substantially enhanced effective mass 
and corresponding small quasiparticlc weight. Generally, the size of the polaron is 
smaller in the non-adiabatic strong-coupling regime due to the faster lattice response. 
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Figure 2. As in figure 1 but for tlie anti-adiabatic regime (7 = 4). 

4.. 1.1. One dimension Figure 1(a) shows results for a single electron in a ID system 
for 7 = 0.1. For K = Q. starting with the non-interacting value ~2t, the absolute 
value of the kinetic energy is reduced with increasing A. The crossover near A = 1 is 
continuous (the ilso applies in higher dimensions [49]). 

Turning on the confinement, we see that results are very similar for small 
K/t = 0.1 (weak confinement). Remarkably, in the case of strong confinement — 
the size of the quantum dot is very small for K/t = 2 — the kinetic energy displays the 
convex behaviour characteristic of the strong-coupling regime already near A = 0. This 
may be understood as a result of squeezing of the polaron state: For weak coupling, 
the free {K — 0) polaron size is larger than the dot size set by so that the kinetic 
energy is determined mainly by the confinement. In contrast, for strong coupling, the 
polaron size is smaller than the dot size and the kinetic energy is almost identical 
to the case K/t = 0. Therefore, all curves tend to the same value — 2t/(l)/7(0) as 
A — > 00. 

The effect of the confinement on the polaron size at the critical coupling A = 1 
(for = 0) is illustrated by the results for the electron-lattice correlation fimction 
Cgp(r) in figure 1(b). For small K, the lattice distortion surrounding the electron 
extends over about two lattice sites (large polaron). The effect of confinement is to 
reduce the polaron size, as is well visible especially for K/t = 2. 

In the inset of figure 1(b) we present data for Cgp(O) as a function of A. As 
pointed out in [50], the small-polaron crossover occurs at the point where the slope 
of Cgp(O) has its maximum, which in the present case confirms the critical coupling 
A = 1 obtained from energetic considerations {K = 0). Furthermore, numerical 
differentiation of the data shown clearly indicates that the critical coupling is reduced 
by about 50 % in the strong confinement regime K/t > 1 owing to the smaller non- 
interacting (A = 0) kinetic energy. 

Figure 2 shows results for the anti-adiabatic regime (7 = 4). For K = 0, the 
kinetic energy in (a) shows a very smooth dependence on the electron-phonon coupling 
strength. The critical value for the small-polaron crossover is given by = 1 (A = 2), 
but no pronounced changes in E^^^ are visible. For K > 0, the form of the kinetic 
energy curve does not change noticeably, because the free polaron size is smaller than 
the confinement length even for g'^ < 1. 

The electron-lattice correlation function in figure 2(b) confirms the smaller size 
of the polaron at the K = critical coupling (cf figure 1(b)). Consequently, the 
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Figure 3. As in figure 1 but for a 2D system. 



dependence on K is much weaker than in the adiabatic regime. As for 7 = 0.1, the 
critical couphng for small-polaron formation is substantially reduced for large K/t, as 
can be inferred from the slope of Cgp(O). 

4.1.2. Higher dimensions For the non-confined case, it is known that the small- 
polaron crossover becomes more abrupt as a function of A in higher dimensions. This 
is well visible from the 2D results for E^^^ shown in figure 3(a), as well as from the 
inset of figure 3(b). 

The effect of confinement is again to reduce the kinetic energy. For K/t = 2, the 
convex strong-coupling behaviour is seen similar to D = 1. An important difference 
to the ID case is that the curves for different K become practically indistinguishable 
for A > 2, whereas in ID (figure 1(a)) such a collapse occurs for A > 3. The origin of 
this effect is revealed in figure 3(b). Obviously, the polaron is smaller in 2D than in 
ID for the same A, reducing the effect of K. The inset again reveals a reduction of 
the critical coupling with increasing confinement. Similar behaviour is also found in 
three dimensions (not shown), with an even stronger dependence of E^^^ and C^p{r) 
on A. 

Let us discuss the relation of our results to previous work. The common conclusion 
is that the confinement demobilises the carrier, which manifests itself either in an 
increase of effective mass, or a decrease of kinetic energy and the polaron size [20-24] . 
Also in accordance with existing work, we find that the polaronic correction to the 
total energy is larger for larger K. For example, AE^ = [E^^^q — -B^=i]/-B^=o ~ ^-^^ 
for K/t = 0.1, whereas AE^ « 1.08 for K/t = 2. 

The reduction of the critical coupling for the existence of a small polaron due to 
confinement complies with the growth of the strong-coupling region at the expense of 
the weak-coupling region in the phase diagram pointed out in [23] . 

Furthermore, we observe a similar dependence of the influence of confinement on 
dimensionality as for the continuum model [22] (in this work the effective mass was 
considered). Comparing the values of the kinetic energy at A = 1, we find for the ratio 
Kini^) / Kini^-^ = 0) = 0-91 {K/t = 0.1), 0.66 {K/t = 1), and 0.53 {K/t = 2) in ID, 
and 0.98, 0.82, and 0.72 in 2D. Hence, confinement effects are significantly stronger in 
ID than in 2D, and this trend extends to D = 3. 
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Figure 4. (a) Kinetic energy Ekin of two electrons in a ID system in the adiabatic 
regime (7 = 0.1) for {7 = as a function of electron-phonon coupling strength A for 
different values of the confinement strength K. (b) Electron-electron correlation 
function C^e (r) as a function of distance r for the same parameters as in (a) and 
A = 0.5. 

4-2. Bipolaron 

Wc now come to the case of two electrons of opposite spin, which can form a 
bound bipolaron state given sufSciently strong electron-phonon interaction. Owing 
to the Coulomb repulsion, the physics becomes much richer, and we again begin our 
discussion with the one-dimcnsional case. 

Compared to the one-electron problem, significantly less work has been done on 
the bipolaron problem in the framework of the Holstein-Hubbard model. Nevertheless, 
from existing studies (see [30,39,40,51] and references therein), the basic physics is 
understood. A summary of (approximate) conditions for the existence of the different 
ground states is given in table 1. 



Table 1. Conditions for the existence of different singlet bipolaron states in 
the one-dimensional Holstein-Hubbard model at weak (VVC) and strong electron- 
phonon coupling (SC) [38,40]. 



U-- 


= 




u >0 




Large 


Small 


Two 


Inter-site 


Small 


bipolaron 


bipolaron 


polarons 


bipolaron 


bipolaron 


A < 0.5 (7 < 1) 


A > 0.5 (7 < 1) 


U > 2Ep (WC) 


U < 2Ep (WC) 




or 


and 






U^2Ep 


g < 0.5 (7 > 1) 


g > 0.5 (7 > 1) 


U > AEp (SC) 


U < AEp (SC) 





4-2.1. One dimension Figure 4 shows results for the simplest case U = 0, in the 
adiabatic regime (7 = 0.1). For K = 0, the two electrons form a bound state if 
AEp > 4t (A > 0.5), where AEp is the bipolaron binding energy in the atomic limit 
and 4:t is the free kinetic energy of the two electrons. The corresponding crossover 
from a large bipolaron (no unbound polarons exist for U = 0) to a small (on-site) 
bipolaron leads to a noticeable decrease of the (absolute) kinetic energy in figure 4(a). 
As for the polaron, the critical value for the crossover may be determined from the 
slope of C (0) or Cgg(O) (see figure 5). 




Figure 5. Electron-electron correlation function Cee(O) as a function of electron- 
phonon coupling strength A for a ID system, 7 = 0.1, and different values of the 
Hubbard repulsion U and the confinement strength K. 



Apart from enforcing polaron effects, confinement enhances the probabiUty for the 
two electrons to occupy the same region of the system (the vicinity of the centre of the 
harmonic potential). This in turn leads to a stronger pairing tendency. Any K/t > 
decreases the kinetic energy in the weak- and intermediate coupling regime. For strong 
confinement, the kinetic energy exhibits the typical strong-coupling dependence on A. 
The reduction of the bipolaron size (the average distance between the two carriers) due 
to confinement is shown in figure 4(b) for the free critical coupling A = 0.5. Finally, 
figure 5 reveals that the critical coupling for small-bipolaron formation is substantially 
smaller in a strongly confined system. 

It is more realistic to consider the case of finite Coulomb repulsion U/t = 4. 
For K = 0, strong-coupling theory predicts a ground state with two unbound (large) 
polarons for 2Ep < U. In the vicinity of 2Ep = U, the so-called inter-site bipolaron 
exists, with the electrons being most likely to occupy neighbouring lattice sites, and 
with only slightly enhanced effective mass as compared to two unbound polarons 
[37,38,40]. The origin of this state are non-Iocal exchange interaction processes [38] 
which allow the particles to gain potential energy from the c;oupling to phonons, and 
at the same time avoid the penalty due to Hubbard repulsion. 

Remarkably, figure 6(a) reveals a critical coupling for the small-bipolaron 
crossover of 2Ep « J7 (A « 1, sec also inset of figure 6(b)), which is in a surprisingly 
good agreement with the anti-adiabatic strong-coupling condition despite 7 = 0.1. Our 
findings extend previous calculations limited either to 7 — 1 [38] or small systems [40] , 
and a phase diagram for the important adiabatic regime will be presented elsewhere. 

As for [/ = 0, confinement leads to a reduction of the kinetic energy and the 
average distance between the two electrons (figure 6(b)). However, owing to the on- 
site repulsion, the results for E^^^^ exhibit weak and strong-coupling behaviour even 
for K/t = 2. 

An interesting point concerns the effect of confinement on the inter-site state. 
The inset in figure 6(b) shows the electron-electron correlation functions C^^(O) and 
Cgg{l). For if = 0, in the vicinity of A = 1, Cgg(l) > C^^{Q) as characteristic for 
the inter-site bipolaron. In contrast, for strong confinement K/t = 2, the region of 
existence of the latter is limited to small A. Clearly, for sufficiently strong confinement, 
the inter-site state will disappear. 

In previous strong-coupling calculations in the framework of the continuum model 
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Figure 6. As in figure 4 but for U/t = A. The inset in (b) shows the electron- 
electron correlation functions Cee(O), Cee(l)- 

[24-26], it has been found that the bound bipolaron state becomes unstable at very 
strong confinement, which was attributed to the increase of the Coulomb interaction 
energy for two spatially close particles. However, this effect has been argued [23, 26] 
to be due to the approximations made, and path-integral calculations [23] as well as 
analytical and QMC calculations [17] suggest that a bound state also exists at weak 
coupling in a confined system. 

Preliminary calculations of the bipolaron binding energy E^{N^ = 2) — 2E^{N^ = 
1) for 7 = 0.1 and U/t = 4 at K/t = 0.1 and K/t = 2.0 in one dimension reveal that a 
weakly bound bipolaron can indeed become unbound due to confinement. This effect 
of on-site Coulomb interaction is expected to be most pronounced in one dimension, 
and this issue will be further investigated in future work. 

It is worth mentioning the relation of our calculations to recent work on molecular 
quantum dots [18, 19] based on models with a single molecular level and a vibrational 
mode, weakly coupled to metallic leads. In such systems, the bipolaron binding energy 
can compensate for the Hubbard repulsion U, giving rise to a net attraction and 
thereby favouring pair tunneling of electrons. Due to the absence of a sizable hopping 
amplitude, the ground-state properties of such a single-site molecule are determined 
mainly by the atomic-limit physics of the present model. In particular, the inter-site 
bipolaron state will be strongly suppressed, and the strong-coupling criteria of table 1 
are expected to hold. The situation may be expected to be different for finite-size 
molecules. 

4-2.2. Higher dimensions Numerical results for the unconfined Holstein- Hubbard 
model in D > 1 have been reported before in 2D for 7 = 1 [39] , and in 3D for classical 
phonons [35], but no reliable results are available for the important adiabatic regime 
7 <C 1. We restrict the discussion to the more general case U > 0. 

As in the polaron problem, the crossover to a small bipolaron becomes more 
pronounced in higher dimensions (note the different ordinate scales in figures 6(a) and 
7(a) and (b)). The critical coupling fov K = from the strong-coupling approximation 
is again set by ?7 ~ 2Ep, corresponding to U/W = X. Indeed, for D = 2 (figure 7(a)), 
we find the crossover near A = 0.5 in agreement with this condition. However, in three 
dimensions, the crossover occurs slightly above A = 0.5, whereas the strong-coupling 
result yields A = 0.33. This deviation may be explained by the larger number of 
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Figure 7. Kinetic energy Ekin for two electrons in the adiabatic regime (7 = 0.1) 
for U/t = 4 as a function of coupling strength A for different values of the 
confinement strength K in (a) two dimensions, (b) three dimensions. Insets show 
the electron-electron correlation function Ceeir = 0) versus A. 



nearest-neighbour sites in 3D, which may enhance the stabihty of extended (inter- 
site) bipolaron states. 

Confinement again leads to a reduction of the critical coupling, as well as to an 
enhancement of bipolaron effects. As for the a single electron, we can calculate the 
change of the kinetic energy with confinement (relative to ii' = 0) at A = 1. This 
yields 0.89 in two dimensions, and 0.95 in three dimensions, so that we can conclude 
that confinement effects are again more noticeable in lower dimensions. 

5. Summary and outlook 

Unbiased quantum Monte Carlo studies of the Holstein-Hubbard model with 

additional harmonic confinement on a discrete lattice have been carried out in one to 
three dimensions, considering the cases of one and two electrons. Technical difficulties 
encountered in simulations using the present world-line method, some of which were 
overlooked in previous work, have been revealed. 

The effect of confinement on the formation of polarons and bipolarons has been 
investigated. Despite considerable simplifications inherent to the model, we believe 
that our exact numerical results arc relevant for quantum dot systems. 

In addition to providing unbiased results for all physically relevant parameter 
regimes, we have for the first time reported on results for a (multi-site) cluster model 
of a quantum dot. The latter is particularly important for the very small dots which 
can be fabricated today, as well as for a description of the small-(bi-)polaron crossover 
associated to a local lattice instability. Quantum phonon effects and electron-electron 
interaction were fully taken into account. 

For one electron, we find that the basic effect of confinement consists in enhancing 
polaronic eff'ects — the polaron size and the critical coupling for the existence of a small- 
polaron both decrease. The influence of confinement is largest in a one-dimensional 
system in the adiabatic regime, and becomes significantly smaller with increasing 
dimension or phonon frequency. In the strongly confined regime, the polaron state is 
squeezed, giving rise to small-polaron physics even for weak or intermediate coupling. 

We have presented the first accurate results for the unconfined Holstein-Hubbard 
bipolaron in the adiabatic regime in more than one dimension. We find indications that 
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confinement can counteract bipolaron formation in the presence of Coulomb repulsion, 

but this issue needs further investigation also in the framework of a model with long- 
range interactions. For finite on-site repulsion, weak-coupling behaviour survives even 
for strong confinement. 

Existing work on (bi-)polarons in quantum dots is almost exclusively restricted 
to equilibrium properties, although most technical applications fall into the non- 
equilibrium regime. Since the fundamental effects of confinement on the ground state 
arc now rather well understood, it is desirable to consider more general situations 
in the future. Apart from effects due to long-range interactions and anharmonic 
confinement in cluster models, this includes many-polaron effects [52], coupling of 
electrons or excitons to light [16], as well as transport properties and time-dependent 
phenomena. Of course, suitable numerical methods will have to be developed to 
address these problems, and work along these lines is in progress. 
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